function [cs, ctab, vmp, voi] = vmp_ClusterTable(hfile, mapno, threshold, opts)
% VMP::ClusterTable  - generate a table with clusters
%
% FORMAT:       [c, t, v, vo] = vmp.ClusterTable(mapno [, thresh [, opts]])
%
% Input fields:
%
%       mapno       map number (1 .. NrOfMaps)
%       thresh      either p/r values (0 .. 1) or t/F value (1 .. Inf)
%                   if not given or 0, uses the LowerThreshold of map
%       opts        optional settings
%        .altmaps   alternative maps to extract values from (default: [])
%        .cclag     flag, interpret the threshold as lag number (false)
%        .clconn    cluster connectivity ('face', {'edge'}, 'vertex')
%        .localmax  break down larger clusters threshold (default: Inf)
%        .localmin  minimum size for sub-clusters (default: 2)
%        .minsize   minimum cluster size (map resolution, default by map)
%        .mni2tal   flag, VMP coords are passed to mni2tal (default: false)
%        .showneg   flag, negative values are considered (default: false)
%        .showpos   flag, positive values are considered (default: true)
%        .sorting   either of {'maxstat'}, 'size', 'x', 'y', 'z'
%        .tdclient  flag, lookup closest talairach label (default false)
%
% Output fields:
%
%       c           1xC struct with properties of clusters
%       t           text table output
%       v           thresholded map (in VMP resolution!)
%       vo          if requested, VOI structure with (TAL) coords
%
% Note: if only one output is requested, the table is text table is
%       returned!!

% Version:  v0.7g
% Build:    9032021
% Date:     Mar-20 2009, 9:38 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 2 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'vmp') || ...
   ~isa(mapno, 'double') || ...
    numel(mapno) ~= 1 || ...
    isinf(mapno) || ...
    isnan(mapno) || ...
    mapno < 1
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end
bc = bvqxfile_getcont(hfile.L);
if mapno > numel(bc.Map)
    error( ...
        'BVQXfile:BadArgument', ...
        'Map number out of bounds.' ...
    );
end

% get bounding box
bb = aft_BoundingBox(hfile);

% get map structure and data
mapno = fix(mapno);
map = bc.Map(mapno);
vmp = double(map.VMPData);

% check threshold argument
if nargin < 3 || ...
   ~isa(threshold, 'double') || ...
    numel(threshold) ~= 1 || ...
    isinf(threshold) || ...
    isnan(threshold)
    threshold = 0;
end
if threshold == 0
    threshold = map.LowerThreshold;
end

% options check
if nargin < 4 || ...
   ~isstruct(opts) || ...
    numel(opts) ~= 1
    opts = struct;
end
if ~isfield(opts, 'altmaps') || ...
   ~isa(opts.altmaps, 'double') || ...
    isempty(opts.altmaps) || ...
    any(isinf(opts.altmaps(:)) | isnan(opts.altmaps(:)) | opts.altmaps(:) < 1 | opts.altmaps(:) > numel(bc.Map))
    opts.altmaps = [];
else
    altmaps = unique(round(opts.altmaps(:)));
    opts.altmapsp = true(1, numel(altmaps));
    opts.altmaps = double(zeros([size(vmp), numel(altmaps)]));
    for amc = 1:numel(altmaps)
        opts.altmaps(:, :, :, amc) = bc.Map(altmaps(amc)).VMPData;
    end
end
if ~isfield(opts, 'cclag') || ...
   ~islogical(opts.cclag) || ...
    numel(opts.cclag) ~=1
    opts.cclag = false;
end
if ~isfield(opts, 'clconn') || ...
   ~ischar(opts.clconn) || ...
    isempty(opts.clconn) || ...
   ~any(strmpci(opts.clconn(:)', {'edge', 'face', 'vertex'}))
    opts.clconn = 'edge';
else
    opts.clconn = lower(opts.clconn(:)');
end
if ~isfield(opts, 'localmax') || ...
   ~isa(opts.localmax, 'double') || ...
    numel(opts.localmax) ~= 1 || ...
    isnan(opts.localmax) || ...
    opts.localmax < 3
    opts.localmax = Inf;
else
    opts.localmax = round(opts.localmax);
end
if ~isfield(opts, 'localmin') || ...
   ~isa(opts.localmin, 'double') || ...
    numel(opts.localmin) ~= 1 || ...
    isinf(opts.localmin) || ...
    isnan(opts.localmin) || ...
    opts.localmin < 1
    opts.localmin = 2;
end
if ~isfield(opts, 'minsize') || ...
   ~isa(opts.minsize, 'double') || ...
    numel(opts.minsize) ~= 1 || ...
    isinf(opts.minsize) || ...
    isnan(opts.minsize) || ...
    opts.minsize < 1 || ...
    opts.minsize > (numel(vmp) / 2)
    if map.EnableClusterCheck
        opts.minsize = map.ClusterSize;
    else
        opts.minsize = 1;
    end
else
    opts.minsize = fix(opts.minsize);
end
if ~isfield(opts, 'mni2tal') || ...
    numel(opts.mni2tal) ~= 1 || ...
   (~isa(opts.mni2tal, 'double') && ...
    ~islogical(opts.mni2tal))
    opts.mni2tal = false;
else
    opts.mni2tal = true && opts.mni2tal;
end
if ~isfield(opts, 'showneg') || ...
   ~islogical(opts.showneg) || ...
    numel(opts.showneg) ~= 1
    opts.negative = (map.ShowPositiveNegativeFlag > 1);
else
    opts.negative = opts.showneg;
end
if ~isfield(opts, 'showpos') || ...
   ~islogical(opts.showpos) || ...
    numel(opts.showpos) ~= 1
    opts.positive = (mod(map.ShowPositiveNegativeFlag, 2) == 1);
else
    opts.positive = opts.showpos;
end
if ~isfield(opts, 'sorting') || ...
   ~ischar(opts.sorting) || ...
    isempty(opts.sorting) || ...
   ~any(strcmpi(opts.sorting, {'maxstat', 'size', 'x', 'y', 'z'}))
    opts.sorting = 'maxstat';
else
    opts.sorting = lower(opts.sorting(:)');
end
if numel(opts.sorting) == 1
    switch (opts.sorting)
        case {'x'}
            opts.sorting = 'z';
        case {'y'}
            opts.sorting = 'x';
        case {'z'}
            opts.sorting = 'y';
    end
end
if ~isfield(opts, 'tdclient') || ...
   ~islogical(opts.tdclient) || ...
    numel(opts.tdclient) ~= 1
    opts.tdclient = false;
end

% some initial checks on threshold
if threshold < 0
    if nargin < 4
        opts.negative = true;
        opts.positive = false;
    end
    threshold = -threshold;
end

% remove invalid entries first
vmp(isinf(vmp) | isnan(vmp)) = 0;

% put some additional fields in opts
opts.mat = bb.QuatB2T;
switch (map.Type)
    case {1}
        thtype = 't-Map';
        thdegf = sprintf('%d', map.DF1);
        opts.tptype = 't';
        opts.tptypedf = map.DF1;
    case {2}
        thtype = 'correlation Map';
        opts.tptype = 'r';
        if bc.FileVersion < 5
            thdegf = sprintf('%d (:= t[%d])', map.DF1, map.DF1 - 2);
            opts.tptypedf = map.DF1 - 2;
        else
            thdegf = sprintf('%d', map.DF1);
            opts.tptypedf = map.DF1;
        end
    case {3}
        thtype = 'CC Map';
        thdegf = sprintf('%d', map.DF1);
        
        % in this case also patch the VMP data, and set fields
        if ~opts.cclag
            vmp = mod(vmp + 2, 1000) - 2;
            opts.tptype = 'r';
            if bc.FileVersion < 5
                thdegf = sprintf('%d (:= t[%d])', map.DF1, map.DF1 - 2);
                opts.tptypedf = map.DF1 - 2;
            else
                thdegf = sprintf('%d', map.DF1);
                opts.tptypedf = map.DF1;
            end
        else
            vmp = floor(0.001 * (vmp + 2));
        end
    case {4}
        thtype = 'F-Map';
        thdegf = sprintf('%d, %d', map.DF1, map.DF1);
        opts.tptype = 'f';
        opts.tptypedf = [map.DF1, map.DF2];
    case {11}
        thtype = 'PSC Map';
        thdegf = 'n/a';
    case {12}
        thtype = 'z(ica) Map';
        thdegf = sprintf('%d', map.DF1);
    case {13}
        thtype = 'Thickness Map';
        thdegf = 'n/a';
    case {16}
        thtype = 'Probability Map';
        thdegf = sprintf('%d', map.DF1);
    case {20}
        thtype = 'Mean Diffusivity Map';
        thdegf = 'n/a';
    case {21}
        thtype = 'Fractional Anisotropy Map';
        thdegf = 'n/a';
    otherwise
        thtype = 'unknown';
end
mthreshp = nan;
if any([1, 2, 4] == map.Type)
    if threshold < 0.2
        mthreshp = threshold;
        switch (opts.tptype)
            case {'f'}
                threshold = custom_finv(1 - threshold, ...
                    opts.tptypedf(1), opts.tptypedf(2));
            case {'r'}
                if opts.negative && ...
                    opts.positive
                    threshold = correlinvtstat(custom_tinv(1 - 0.5 * threshold, ...
                        opts.tptypedf), opts.tptypedf + 2);
                else
                    threshold = correlinvtstat(custom_tinv(1 - threshold, ...
                        opts.tptypedf), opts.tptypedf + 2);
                end
            case {'t'}
                if opts.negative && ...
                    opts.positive
                    threshold = custom_tinv(1 - 0.5 * threshold, opts.tptypedf);
                else
                    threshold = custom_tinv(1 - threshold, opts.tptypedf);
                end
        end
    else
        switch (opts.tptype)
            case {'f'}
                mthreshp = 1 - custom_fcdf(threshold, ...
                    opts.tptypedf(1), opts.tptypedf(2));
            case {'r'}
                if opts.negative && ...
                    opts.positive
                    mthreshp = correlpvalue(threshold, opts.tptypedf + 2);
                else
                    mthreshp = 0.5 * correlpvalue(threshold, opts.tptypedf + 2);
                end
            case {'t'}
                if opts.negative && ...
                    opts.positive
                    mthreshp = 2 - 2 * custom_tcdf(threshold, opts.tptypedf);
                else
                    mthreshp = 1 - custom_tcdf(threshold, opts.tptypedf);
                end
        end
    end
end
thead = cell(1, 5);
thead{1} = sprintf('   Clustertable of map:   "%s"', map.Name);
thead{2} = sprintf('           Type of map:   %s', thtype);
thead{3} = sprintf('    Degrees of freedom:   %s', thdegf);
thead{4} = sprintf('   Cluster k-threshold:   %d mm^3 (%d voxel)', ...
    opts.minsize * prod(bb.ResXYZ), opts.minsize);
thead{5} = sprintf(' Applied map threshold:   %.5f (p < %.5f)', ...
    threshold, mthreshp);

% pass on to clustervol
[ctab, cs, vmp] = clustervol(vmp, threshold, opts.minsize, opts);

% store back in struct if cluster-threshold enabled
if opts.minsize > 1
    bc.Map(mapno).VMPDataCT = dilate3d(vmp ~= 0);
    bvqxfile_setcont(hfile.L, bc);
end

% extend ctab
ctab = sprintf('%s\n%s\n%s\n%s\n%s\n\n%s', thead{:}, ctab);

% create VOI structure ?
if nargout > 3
    voi = bless(BVQXfile('new:voi'), 1);
    voic = bvqxfile_getcont(voi.L);
    if ~isempty(cs)
        voic.VOI(numel(cs)).Name = '';
    else
        voic.VOI(:) = [];
    end
    for cc = 1:numel(cs)
        voic.VOI(cc).Name = sprintf('Cluster%04d_%d_%d_%d_%s', ...
            cc, cs(cc).rwcoords(1, :), ...
            strrep(regexprep(char(cs(cc).talout), '\(.*$', ''), ' ', '_'));
        if cs(cc).values(1) >= 0
            voic.VOI(cc).Color = floor([255, 127.999 * rand(1, 2)]);
        else
            voic.VOI(cc).Color = floor([127.999 * rand(1, 2), 255]);
        end
        voic.VOI(cc).NrOfVoxels = size(cs(cc).rwcoords, 1);
        voic.VOI(cc).Voxels = cs(cc).rwcoords;
    end
    voic.NrOfVOIs = numel(voic.VOI);
    bvqxfile_setcont(voi.L, voic);
end

% return table if only one output is requested
if nargout < 2
    cs = ctab;
end
